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ABSTRACT 

In this paper we introduce the concept of Direct Statistical Simulation (DSS) for 
astrophysical flows. This technique may be appropriate for problems in astrophysical 
fluids where the instantaneous dynamics of the flows are of secondary importance to 
their statistical properties. We give examples of such problems including mixing and 
transport in planets, stars and disks. The method is described for a general set of 
evolution equations, before we consider the specific case of a spectral method optimised 
for problems on a spherical surface. The method is illustrated for the simplest non- 
trivial example of hydrodynamics and MHD on a rotating spherical surface. We then 
discuss possible extensions of the method both in terms of computational methods and 
the range of astrophysical problems that are of interest. 



1. Introduction 



The modeling of astrophysical phenomena is often limited by the huge range of spatial and 
temporal scales that need to be resolved in order to describe accurately the dynamics. In many cases 
the large-scale behavior of cosmic bodies depends on interactions at smaller scales that need to be 
represented properly for a complete understanding of the astrophysical phenomenon in question. 
The situation is usually complicated by the requirement of including the back-reaction of the 
large-scale environment on the smaller scale dynamics in a self-consistent manner. These types of 
problems are ubiquitous in astrophysics, but we list here some important examples. The transport 
of angular momentum in accretion disks may be mediated by the (magneto-)rotational turbulence 
that is present in the disk (see e.g. Balbus & Hawley 1998). This turbulence is itself driven and 
modulated by the large-scale environment of Keplerian rotation and large-scale magnetic fields (see 
e.g. Jamroz et al.|2008 ). The differential rotation pattern in stars (including the sun) arises through 
an interaction of buoyancy-driven turbulence and rotation, with Reynolds stresses at intermediate 
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scales leading to correlations that drive large-scale flows that themselves act back on the turbulence 



(Ruediger 1989 Brun & Toomre 2002 Rempel 2005 Miesch 2005). This situation is mirrored in 



planets, where convective processes may create stresses leading to large-scale flows. Such stresses 
create turbulence in stably stratified outer weather layers (for example in Jupiter and Saturn) that 
may drive the formation of jets (see e.g. Scott fc Polvani||2008 and the references therein). 



There are a number of approaches to modeling the fluid interactions in astrophysical objects. 
The approach taken often depends on whether it is the dynamics or statistics of the system that is of 
interest. Sometimes information about the dynamics — that is the precise evolution of a particular 
realisation of a system — is sometimes required for prediction or to compare with observations. 
More likely is that the statistics — i.e. the average properties of an ensemble of evolutions — is of 
interest; this may give more insight into the underlying physics of the system. 

Theoretically and computationally, a natural first approach is to perform direct numerical 
simulations (DNS) of the fluid (or MHD) equations for the system. This approach is the most 
straight-forward and has led to breakthroughs in many branches of astrophysical fluid dynamics. 
This approach lends itself naturally to determining the dynamics of a given system. However the 
extreme nature of the astrophysical turbulent environment ensures that not all spatial scales may 
be faithfully represented even on the most massive parallel computers available today. For this 
reason the practitioners of DNS must accept that they are not in the correct parameter regime 
or may claim that the parameters take into account the effects of scales below the grid cut-off 
via eddy diffusivities (sometimes termed turbulent transport coefficients). These diffusivities are 
usually chosen in a plausible but ad-hoc manner. Moreover, DNS may not be an efficient algorithm 
for determining statistics, since the ensemble over a large number of expensive calculations may be 
required in order to achieve meaningful statistics. 

An alternative approach, which is not useful for determining dynamics but may be useful 
for statistics, is to derive evolution equations for the large-scale dynamics and to formulate closure 
models for net effects of the dynamics at moderate and small scales. Such models have a long history 



in astrophysics and have also achieved some measure of success (see e.g. Kitchatinov & Ruediger 



(1995 


); Ogilvie 


(2003); Rempel 


(2005)). This approach often utilises (either implicitly or explicitly) 


moment hierarchies (see e.g. ( 


]anuto et al. (1994 


2001 


); 


Farrell &; Ioannou ( 


2008): 


Garaud et al 



(2010)). In particular it is customary to relate the average of local interactions of the small scale to 
the local values of the large-scale fields. The weakness of such models is that they usually rely on 
some ad-hoc assumption to close the model — parameterising the interactions between large and 
small scales — and often make the assumption of homogeneity or isotropy. They sometimes find it 
difficult to include self-consistently the effects of the dynamic large-scale environment. Often it is 
the case in astrophysics that regions of strong transport lie in close proximity to regions of weak or 
no transport or mixing — for a nice example see the jets in Jupiter — and so closures that rely on 
homogeneity may lead to misleading large-scale dynamics. Moreover it is often the case that the 
inclusion of such closure models introduces new adjustable parameters to the problem that can be 
tuned to fit observations and that little is known about the sensitivity of the large-scale dynamics 
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to changes in the parameterizations of the scales not captured. 

In this paper, we present a new approach to the problem of describing astrophysical flows with 
a range of spatial scales, which we believe will prove useful for a certain class of problems with 
large-scale inhomogeneous and anisotropic flows. Specifically we describe the development of effi- 
cient numerical algorithms to solve truncated hierarchies of cumulant equations, leading directly to 
the statistical description of astrophysical flows, and we show that this direct statistical simulation 
(DSS) is able to reproduce qualitatively statistics obtained by time averaging 

DNS^ That DSS 

has several advantages over DNS has long been recognized, going back at least as far as a seminal 



monograph by Lorenz (1967). Low-order statistics are smoother in space and stiffer in time than 



the underlying detailed flows. Statistically stationary fixed points or slowly varying statistics can 
therefore be described with fewer degrees of freedom and also can be accessed more rapidly. Con- 
vergence with increasing resolution can be demonstrated, obviating the need for separate closure 
models of the subgrid physics, although these may be included in a natural statistical framework. 
Finally and most importantly, DSS leads more directly to insights, by integrating out fast modes, 



leaving only the slow modes that contain the physical information of most interest (Lorenz 1967 



Marston 2010). 



In this paper we develop the techniques illustrated recently in Marston et al. (20081, where 



the prototypical problem of barotropic flow relaxing toward a point jet was considered and the 
statistics obtained by DNS were found to be in good qualitative agreement with those found from 
a second-order cumulant expansion. We begin by examining the general case of constructing a 
cumulant hierarchy for the evolution of a number of dynamic variables. We describe the derivation 
and solution of the cumulant equations for the general case, before focussing the discussion on the 
case of spherical symmetry, where computational efficiencies are available. 

Having described the method in general, we illustrate the advantages of the method for a simple 
model of the interaction of turbulence and mean flows that may be relevant to the generation of 
zonal flows in stable layers in planets and stars. The model describes the two-dimensional evolution 
of flows and magnetic fields on a spherical surface. Such an evolution is non-trivial as it is known 
that for the hydrodynamic problem the Reynolds stresses act to drive inhomogeneous zonal flows; 
this type of behavior is difficult to parameterize in sub-grid scale closures. These models and their 
generalizations have been used to describe the dynamics of the outer layers of giant planets such 



as Jupiter (see e.g. Scott & Polvani (2008) and the references therein) though competing theories 



for the generation of zonal flows via deep-seated convection (see e.g. Jones & Kuzanyan (2009) and 



the references therein) are also available. Furthermore there has been much interest in the MHD 
version of this problem owing to its importance in the dynamics of the solar tachocline (see e.g. 



Tobias et al. (2007); D. W. Hughes, R. Rosner, fe N. O. Weiss (2007)) and potentially in the outer 



layers of extra-solar planets (Staehling Sz Cho (2006)). Both of these environments are believed to 



1 We note here that we choose the terminology direct statistical simulation as we are solving for the statistics 
directly. 
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be turbulent, stably stratified and magnetized. 



The tachocline is believed to play a crucial role in the generation of the eleven year solar cycle 



(see e.g. Tobias & Weiss (2007)) and angular momentum transport through the tachocline may be 



responsible for spinning down the solar interior. One crucial issue to be resolved is therefore the role 
of turbulence in transporting angular momentum in the tachocline, and this has been addressed 



both in the hydrodynamic and magnetohydrodynamic settings (see e.g. Spiegel & Zahn (1992); 



Gough & Mclntyre| fll998| ); |McIntyre| (|2003j); |Tobias et al.| (|2007j)). What has been shown is that 
whereas anisotropic hydrodynamic two-dimensional turbulence leads to the efficient formation of 
zonal flows via Reynolds stresses; the addition of a magnetic field leads to Maxwell stresses that can 
oppose the formation of jets. The suppression of jets is a function of the strength of the large-scale 
magnetic field and the local magnetic Reynolds number Rrn- 

The paper is organised in the following manner: in the next section we introduce the general 
method and the computational savings that can be achieved for the case of spherical symmetry 
in two dimensions. In section 3 the particular model of MHD turbulence on a rotating spherical 
surface is introduced and a comparison of the large-scale dynamics of DNS and DSS is mad^J We 
conclude by discussing extensions to the method and speculating on the range of problems where 
such a technique may be of use. 



2. Formulation of the Model 

In this section we describe the derivation of a general fully spectral algorithm for the direct 
statistical simulation of astrophysical flows [^} We develop the method for the typical case of 
equations with quadratic nonlinearities, before specialising to systems with spherical symmetry in 
two dimensions. 

Consider a system that is represented by partial differential evolution equations (PDEs) for a 
number r of scalar fields. Typically such a system may be solved directly by discretising the PDEs 
using a finite-difference, finite volume or finite element method or by deriving equations for the 
amplitude of modes in a spectral expansion. Formally, this transforms the PDEs into a finite set of 
ordinary differential equations (ODEs) that may be integrated forward in time. If the discretisation 
is performed at s discrete points (or for s spectral modes) then the evolution equations can take 



2 A more in-depth discussion of the dynamics, including the calculation of turbulent transport coefficients is 
included in a forthcoming paper. 

3 In a future paper we shall describe the adaptation of pseudospectral methods for DSS 
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the form 

en = Ai+ Bij qj + C ijk qjq k + /j(i) 
</*(*)> =0 

(mf J (t'))=r t ,5(t-t'). (i) 

where 1 < i, j, k < rs and Ai, and C^ k are the coefficients. Here the qi are the discretised 
values of the dependent variables (or the amplitudes of the relevant spectral modes) ; typically these 
represent a vector of the values of the fluid properties. We also note here that there is an implicit 
sum over repeated indices. 

Hereinafter, to fix ideas, we shall think of the qi as representing the amplitudes of the spectral 
modes of a vector of dependent variables — and shall give a concrete example in the next section. 
The forcing fi{t) can then be interpreted as the statistical forcing of the relevant spectral mode. 



2.1. The general case 

2.1.1. Reynolds decompositions, cumulants and moments 

One way to formulate the cumulant expansion is by carrying out a Reynolds decomposition of 
the dynamical variable qi into the sum of a mean value and a fluctuation (or eddy): 

Qi = {qi)+q'i with(^) = (2) 

where we defer, for now, choosing the type of averaging denoted by the angular brackets (). Typical 
choices are temporal or zonal averages]^] or averages over an ensemble of initial conditions or an 
ensemble of realisations. 

Once the Reynolds decomposition has been implemented, progress is made by defining the first 
three equal-time cumulants q, Cij, and c^k of the combined scalar fields (%) as[^j 

a = (Qi) = m i, 
Cij = {qWj)i 

= (QiQj) ~ (Qi)(qj) = m ij - m i m h 

Cijk = (q'i,q'jq'k)i 

= (qiqjqk) - {(qi)(qjqk) + perms) + 2(q i )(q j )(q k ), 

= m ijk - m,im jk - mjm ik - m k m,ij + 2m;mjm fc , (3) 



4 Although the qt are functions of time only, a zonal average may be taken by keeping only those modes that 
correspond to axisymmetric modes (for an expansion in spherical basis functions this is equivalent to keeping the 
m = modes). 

5 These definitions are sufficient for cumulant hierarchies truncated at either second or third order; for higher order 
hierarchies corresponding definitions of the higher order cumulants are required. 



- 6 - 



where mi, rriij, and rriijk are respectively the traditional definitions of the first, second and third 
moments. We stress here that the second and higher cumulants contain information about corre- 
lations that are non-local in space and therefore include interactions that are not included in the 
simple local moment hierarchies discussed in the introduction. For this reason this approach is 
more tailored to inhomogeneous problems. 



2.1.2. Derivation of the cumulant hierarchy: the H op f functional approach 

The hierarchy of equations of motions for the evolution of the cumulants can be obtained 
directly be differentiating Eqs. [3] with respect to time and using Eqs. [TJ together with repeated 
back substitution. A more elegant method is to introduce variables pt that are, in analogy to 
quantum mechanics, conjugate to the qi in the sense that % = —id/dpi as in Eq. [7] below ( pafc 
Marston 2005). Then one may define the Hopf generating functional ( Frisch|[l995[ ): 



p] = e«), (4) 

recalling the summation over repeated indices. The Hopf functional obeys a Schrodinger-like equa- 
tion: 

zJU = H% (5) 

with linear operator H given by: 

/ d d 2 \ 

H = pi I —Ai + %Bij — + C ijk , (6) 

V dpj opjdpkj 

as can be verified by combining Eqs. |4j[5j and [6] to reproduce Eq. [I] in the absence of any stochastic 
forcing. 



As Eq. |5| is linear in \P, the average ^ = fi?[q(t), p]) obeys the same equation; however ^ 
encapsulates information about the equal-time moments, as can be seen by repeated differentiation 
of Eq. [4] with respect to pi, followed by averaging: 



dp h dp i2 ■ ■ ■ dp in 



(7) 

Pi=0 



(For the case of time-averaging, the statistics do not vary in time, = 0, and the statistics are 
obtained from the solution of the time-independent equation = 0.) The Hopf functional ^ 
may also expressed as the exponential of a power series in pi, the coefficients being the cumulants: 

^ = exp I ia(t) pi - 2jCy(£) PiPj ~ ^<Hjk{t) PiPjPk + ■ ■ ■ \ (8) 
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as can be checked by use of Eq. [7] to reproduce the moments in terms of the cumulants, Eqs. |3j 
Stochastic forcing can now be included with the addition of the T{j term: 



H 



Pi 



-Ai + iBi 



_d_ 



+ C, 



2 



Opjdpk 



(9) 



Upon substituting Eq. [8] into Eq. [5] and collecting powers of pi one obtains the equations of motion 
(EOM) for the cumulants that truncated at third order read 



Cj — Ai + Bij Cj -\- Cijk {pj C-k pj'fe) 

tij = {2Bik Ckj + Cike (4q Cjk + 2cjke)} + Tjj 

Cijk = {3Ba C£jk + GCk£ m {Cijm Q + Cj m Cjt)} — /i Cijk + 0(Cijkl) 



(10) 



where we defer discussion of the parameter [i until later. Here for compactness we have introduced 
the short-hand notation {} to denote symmetrization over indices 



{2B ik c^} = Bik Ckj + Bjk Cki 
that maintains symmetries Cy = tji and similarly for the third cumulant. 



Truncated at second-order (CE2) the cumulant expansion is realizable (Salmon 



1998 



(11) 



and well- 



behaved in the sense that the energy density is positive and the second cumulant obeys positivity 
constraints. Going to third-order (CE3) and beyond introduces difficulties. A phenomenological 



eddy-damping parameter (Orszag 1977 Andre 1974) \i that models the neglect of the fourth and 



higher cumulants from the hierarchy is included in the last of Eq. 10 and is required to prevent 
blow-up. This ad-hoc procedure is somewhat unsatisfactory and more robust methods may be 
necessary. Indeed determining reliable methods of truncating the hierarchy is a matter of current 
research. 



2.2. Symmetry and the derivation of reduced cumulant equations. 

In principle, the general set of cumulant equations in Eq. [TU] can be solved with enough com- 
putational effort. However, efficient algorithms can be developed if the underlying system exhibits 
further symmetries. This is typically the case for astrophysical systems, which usually exhibit 
spherical or cylindrical symmetry or a corresponding translational symmetry in a local Cartesian 
domain. We discuss in detail here the case of cumulants in a sphere. 



6 CE2 can be viewed as the exact solution of a stochastically-driven linear model 
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2.2.1. Equations for fully spectral DNS on a sphere 

For systems with an underlying spherical symmetry, the spectral expansion of the dependent 
variables discussed in section [2] often takes the form 

L;M 

q = ^q im (r)Yr(0,(f>), 

£;m 

= EqUr) (-ir/ ^ 1 ^^; PT^ey^, (12) 

where r is spherical radius, 6 is co-latitude and (j) is longitude. Here the qi m (r) are complex 
functions and the PY 1 are associated Legendre functions. Furthermore on a spherical surface the 
r-dependence is absent and a fully spectral representation of the equation of motion (equation [I]) 
can be written as 

qi m = Ag 5 m fi + Bt + ftm(t) 

ii 

m=m\+rri2 



m=m\ —m2 



£l,t2,mi,m 2 

We note here that, because the scalar fields are real- valued in coordinate space, we may focus on 
the evolution of modes with m > as modes with m < may be obtained by complex conjugation. 
Moreover, for simplicity in the above and in subsequent equations the index that encodes which 
state variable is being solved for has been subsumed into the £ label. 

The quadratic nonlinearities have their origin in the Jacobians of Eqs. 16 with coefficients 
representing amplitudes for the scattering of two waves with m > 0; C (_) are for waves with 
m > and m < to scatter. The amplitudes of these coefficients are constructed from the matrix 
elements of the Jacobian: 

lSL-^2 = / d0 PPioosO) P^(cos8) ^P^(cos9) (14) 

where m = mi ± m.2. Integrals 1^ are evaluated in a numerically exact manner by Gaussian 
quadrature. 



2.2.2. Equations for fully spectral DSS on a sphere 

Similar considerations lead to an efficient representation of the cumulant hierarchy for a spher- 
ical shell. These considerations can then be combined with the knowledge of the underlying sym- 
metries of the statistics themselves to derive reduced hierarchies of cumulant equations. These 
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symmetries are preserved whether zonal, temporal or ensemble averages are used. Statistics on 
the rotating sphere exhibit azimuthal symmetry. The simplest conceptual choice for the averaging 
operation () is therefore the zonal average and we choose that here, and then follow that with a 
running time average. On symmetry grounds, the first cumulant must be independent of longitude 
4> and therefore in the spherical harmonic basis only the m = mode q = {qe, m =o) is non-zero. 
Similar symmetry arguments yield the result that the second cumulant depends on the latitudes of 
the two field points, but only on the difference between their longitudes. It can therefore be written 
as Q 1 £ 2m = {qe im Q£2-m) — ci 1 ce 2 S m o- Furthermore, zonal averaging then requires that Q 1 ^ 2m= o = 0. 
Similarly the third cumulant is a function of only 5, not 6, wavenumbers, i.e. it can be written 
as ce ; i imVi i 2m2 . Moreover, because the scalar fields are real- valued in coordinate space, we have 
ce 1 e 2 m = c e 2 £ im - F° r models with an imposed north-south reflection symmetry about the equatoiQ 
the cumulants respect further constraints: q vanishes for all even £ and Q 1 £ 2m = if l\ is odd and 
£2 is even, and vice-versa. All of these symmetries therefore lead to a computational saving. 

We consider here the simplest non-trivial case where the hierarchy is truncated at second order 
(CE2), i.e. all higher cumulants are set to zero. The EOM for the cumulants in the basis of spherical 
harmonics are then 

+ B h 

+ Cetera C ? C «2»n + C e£eO;i'rn C( - C hi'm, (15) 

where again the convention of summation over repeated indices has been adopted. (There would 
also be a contribution to the first cumulant from C^] -^'0 C( - Cp ^ vanishes for the problems 
considered here as the Jacobian of two fields with no longitudinal dependence is zero.) These 
equations may be compared to a coordinate-independent version given by Equations (21) and (22) 
in Marston et al. (2008). That only the eddy- mean flow interaction is retained in CE2 may be 



seen by noting that the coupling of the first cumulant with the second involves no mixing of the 



azimuthal wavenumber m (only a single m appears in Eqs. 15). Eddy-eddy scattering occurs only 
at third and higher orders^! 



3. Turbulent driven MHD on a Spherical Surface: The Model 



We consider a simple two-dimensional model of a stably stratified region of hydrodynamic 
or MHD turbulence. This is the simplest extension of the local /3-plane model considered by 



Tobias et al. (2007). We stress again that, although this system is of interest in its own right and 



the interaction of Reynolds and Maxwell stresses play an important role in the dynamics of the 



7 this is not necessary, but is computationally expedient 

8 We shall investigate including higher orders in the hierarchy for the cumulant expansion in subsequent papers. 
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tachocline and other regimes of stably stratified MHD turbulence, in this paper we are utilising 
this model as a non-trivial example of the utility of direct statistical simulation. We therefore defer 
discussion of the interaction of the stresses for a subsequent paper. 

The behaviour of such a system in two dimensions can be described by the evolution of two 
scalar fields, namely the relative vorticity ((9,(f>,t) (with 9 being co-latitude and cfi longitude, as 



before) and the scalar potential for the magnetic field A(9,(p,t) (cf Tobias et al. 2007). When 
extended to the sphere rotating at angular rate f2 these may be written: 

q = J[q, 4>} + J[A,V 2 A'}- K (-v 2 V\ + f(t), 

A = J[A, rp] +rjV 2 A', (16) 

where on the unit sphere the Jacobian is given by 

1 fdqdip dqdi/j\ 

J[q ' ® = M\Mde~d6W' (17) 

Here 

C = v 2 ^, 

q = C + 2ftcos0, 

A = A' + B cos9. (18) 

Hence q is the absolute vorticity. Here k is a frictional term, u 2 is a hyperviscosity whilst f(t) is the 
stochastic forcing. Magnetic diffusion is explicitly included through a magnetic diffusivity rj, since 
we believe it is important to capture this process correctly. We note here that the parameter Bq 
measures the strength of a toroidal imposed magnetic field, which is held fixed in time and note that 
such a field can not be self-consistently maintained in a strictly two-dimensional calculation. The 
equations have been scaled so that the magnetic field is measured in units of the Alfven velocity. 
For purely hydrodynamic simulations we simply set Bq = A = 0. 

These equations may be then be written in the form of equation (with r = 2) by setting 
the absolute vorticity and magnetic potential scalar fields into two layers labelled by q a with qi = q 
and q 2 = A and discretising the system either to obtain equations for the spectral amplitudes of 



the form equation (13) or more conveniently for computation a finite difference representation on a 



spherical geodesic grid. In a similar manner the spectral representation of the cumulant equations 



(i.e. equations ( 15 )) can simply be derived once the coefficients in equation ( 13 ) have been calculated 
for this model. 
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4. Comparison of DNS and DSS 
4.1. Numerical implementation of DNS and DSS 

4-1.1. DNS 

Direct numerical simulation of the two-dimensional system has been implemented using two 



different techniques. The EOM given by equation (13) may be integrated forward in time in their 
pure spectral form using a standard fourth-order accurate Runge-Kutta algorithm with an adaptive 
time step, though in practice it is much faster work directly in real space on a spherical geodesic 
grid as we do here. The fully spectral code is therefore only used as a validation of the geodesic 
code below and a useful comparison with the fully spectral direct statistical simulation. 

The most efficient numerical integration of the DNS EOM is carried out in real space on a 
spherical geodesic grid ( Heikes fe Randall|1995 ) of D cells with the use of the second-order accurate 



leapfrog algorithm and a Robert filter. A multigrid algorithm solves Poisson's equation at each 
time step. 



4.1.2. DSS 



We take advantage of the stiff nature of the spectral EOM for the cumulants (equations 15) 



These are integrated forward in time using a semi-implicit backward Euler Full Orthogonal Method 



(Saad 2003) that is based upon Krylov subspaces and that permits a much longer time step than 



is possible for explicit integration methods. 



We note here that integration of the EOM for CE2, equations (15), requires of order L S M 
operations at each time step, where < I < L and < m < min(^, M) define the spectral cutoffs. 
A pseudo-spectral implementation of the EOM would require the same order of operations on the 
sphere and thus offers no advantages over the pure spectral method used here. We find that all 
Qi^ 2 m with m greater than the maximum azimuthal wavevector of the stochastic forcing vanish, 
hence the spectral expansion can be severely truncated by restricting M <C L without loss of 
accuracy. This results in substantial speed-up and a reduction in the required memory. Moreover, 
only a subset of the possible coefficients of the quadratic nonlinearities, Cj.^ m .^ aTn and Cj^jo^'m' 



with 4 indices appear in Eq. 15 resulting in reduced memory usage. 



Finally we note that the code implementing both DNS and DSS (via CE2) is written in the 
Objective- C++ programming language and runs on Apple computers (OS X 10.6) utilizing C- 
blocks and grand central dispatch (gcd) for efficient SMP parallelism. We stress that the DSS can 
run an order of magnitude or more faster than DNS. 
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4.2. Conservation Laws, Model Parameters and Initial Conditions 



In the absence of damping and driving forces, the EOM for the cumulants, like the EOM for 
the vorticity and magnetic potential have a number of conservation laws. For example, in the 
hydrodynamic case, kinetic energy, enstrophy and angular momentum are conserved, whilst for the 
MHD case the conserved quantities are angular momentum, total energy, cross-helicity, and the 
mean squared potential. Moreover, for stochastic forcing restricted to wavevectors \£\ > 0, the case 
considered here, the angular momentum in the CE2 remains exactly zero, in contrast to DNS. 

Just as for direct numerical simulations utilising spherical harmonics there are convenient 
expressions of the average values of various quantities in terms of the low-order cumulants. For 
example the mean cross-helicity is given by: 

= -„ y^jcmem + cu c 2 e S m0 ) (19) 

tm 

where the two layers are labelled explicitly in the final line. Similar expressions are available for 
the averages of other quadratic quantities. 

The models are formulated on the unit sphere with a timescale such that the sphere complete 
a full rotation in one day of model time. All model parameters may be defined in terms of these 
length and time scales; for instance Q = 2tt. Friction removes energy at long length scales and 



is parameterized by rate k. The hyperviscosity v 2 that appears in Eq. (16) is included solely to 
absorb enstrophy at the smallest resolved scales. Consequently it is rescaled with the grid size or 
spectral cutoff so that 

v 2 — > f2*(D/4)~ 2 geodesic, 

v 2 -)• v 2 * (L(L + l)y 2 spectral, (20) 

where the maximum eigenvalue of —V 2 on a geodesic grid with D cells is approximately D/2. Thus 
for u 2 = I (the case we consider here) features on the smallest length scales are dissipated on a 
time scale of order 1 day. 



Stochastic forcing is confined to wavevectors L m i n < t < L max and M m i n < \m\ < L 



max ■ 



Within this range of wavevectors the forcing fi m that appears in Eq. ( 13 ) is given by 

f im {t) = v f F/A * gaussian(t/A) (21) 

where gaussian(f) is a complex number randomly drawn, for each value of I and m, from a normal 
distribution of zero mean and unit variance that smoothly transitions from one random number to 
the next over a time period of A. We set A = 0.1 which is large compared with the time step, but 



small compared with advective time scales. Consequently in Eq. (15) we have 



tm 



fov L mm 

<£< L max and M min < \m\ < £ 

{22) 

else. 



- 13 - 



In the following we hold fixed k = 0.02, V2 = 1, F = 0.2, and (for the magnetic cases) 77 = 10 4 . 
We study the evolution of the systems (DNS and DSS) for two different choices for the range of 
the forcing wavevectors {L min ,M min ,L max }. 

We close our description of the set-up of the models by commenting on the choice of initial 
condition. The DNS integrations are started from rest with zero perturbation to the imposed field. 
For the DSS, at the start of the CE2 integration we set the first cumulant q = and the second 
cumulant c^ 2m = C2 <5^ 2 which corresponds to initial short-ranged correlations in the vorticity. 
At low resolutions, the fixed point sometimes has jets that move in directions opposite to those 
found in DNS; this fixed point, which is an artifact of the spectral truncation, can be avoided by 
initializing q with small values. 

5. Results 

5.1. Small-scale forcing: L m i n = 8, M m i n = 8, and L max = 12 

We begin by considering the hydrodynamic and magnetohydrodynamic evolutions for the case 
where the system is forced solely at small scales in the vorticity equation. The DNS of the hydro- 
dynamic case is performed until a statistically steady state is reached and meaningful statistics can 
be calculated. For this case, this has occurred by t ~ 1000; after this time a running time-average 
is performed for another 1000 days. Here the small-scale driving leads to the formation of flows 
on a range of scales including large-scale jets as shown in Figure [I] (top panels) which show the 
instantaneous relative vorticity (left) and relative zonal velocity (right). These clearly show the 
formation of a prograde (westerly) jet at the equator with two retrograde (easterly) jets at high 
latitudes, with the total angular momentum of the fluid remaining close to zero. As we shall see, 
these jets are driven by the flows on smaller scales. 

The history and statistics of these hydrodynamic jets for DNS is displayed in the timelines in 
the upper panels of Figure [2j In the left portion of each panel the relative vorticity and relative 
zonal velocity (averaged over a period of 10 days) is shown as a function of latitude and time. At 
1000 days (half-way through the evolution — signified by a vertical line in the figures) temporal 
averaging is switched on and a running average from that point is displayed in the figures. This 
running average eventually settles down to show the mean position and strength of the jets. 

Figure [2] compares these timelines with those calculated by DSS for the same parameter values. 
The direct statistical simulation achieves remarkable agreement with the DNS in both the position 
and strength of the jets. This is confirmed in Figure [3] which demonstrates that the time-averaged 
zonal mean zonal velocity as calculated by DSS agrees well with DNS except at high latitudes. 
Moreover, whilst the DSS respects the north-south symmetry as expected, for the DNS the average 
position of the prograde jet is slightly off-equator, reflecting the finite length of data over which 
the averages are calculated. Figure [3] also demonstrates that good convergence with increasing 
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resolution is achieved, both for DNS and for DSS. 



That DSS and DNS agree well is reflected in the data recorded in Table 1. There we give some 
non-dimensional ratios that can only be calculated once the kinetic and magnetic energy are in a 
statistically steady statej^] These are defined as 



Ro 



(u 2 + v 2 ) 2 



Rm 



(u 2 + v 2 )zL 



(23) 



where we recall that $7 = 2ir, L = 1 and r\ = 10 4 . That DSS is able to reproduce the jets using 
a cumulant hierarchy truncated at second order is an interesting result. It is evidence that the 



forward enstrophy cascade and anisotropic backward energy cascade (Kraichnan & Montgomery 



1980 , Salmon|1998 ), which is frequently invoked to explain their existence, is in fact not necessary — 
there can be no cascade in the absence of eddy-eddy interactions. We note here that the enstrophy 



cascade argument has also been questioned in the context of planetary atmospheres (Vallis 1992 



Schneider Sz Walker 2006 O' Gorman Sz Schneider 2007) where shearing and modification of the 



thermal structure of the atmosphere by eddy fluxes weaken the eddy-eddy interactions. Here it is 
therefore Reynolds stresses that are primarily responsible for the build-up of the zonal flows. This 
result is important for our understanding of the driving of zonal flows in planetary atmospheres. 

We now examine the effects of including a toroidal magnetic field and examine the dynamics 
for two imposed field strengths, Bo. For B$ = 0.1 the onset of jet formation is delayed but in both 
DNS and DSS the system eventually settles into a statistically steady state as shown in Figure |4j 
For both methods, for this relatively weak imposed mean field, there is eventually little suppression 
of the jets, with slightly more suppression occurring in the DSS. For this choice of parameters, the 
magnetic energy is small compared with the kinetic energy of the flow (4% for both DNS and DSS 
(CE2)), and so it is to be expected that the role of the magnetic field will be secondary. Moreover the 
magnetic field has been expelled to high latitudes by the strong jets and turbulence at low latitudes 
(see Figure [7]). This flux expulsion (Weiss 1966 Tao et al. 1998) leads to separated regions with 
different dynamics; at low latitudes, where the field is weak, the hydrodynamic evolution continues 
unimpeded, whilst at high latitudes the magnetic field leads to some suppression of the jets. 

At Bq = 0.5, however, strong qualitative changes are plainly evident as the jets are destroyed 
by the fluctuations in the magnetic field, both in DNS and in DSS, in agreement with the findings 



of Tobias et al. (2007). Small remnants of the jets persist in DSS at high latitudes, where the 



imposed toroidal field is weakest (see Fig. [5]). DNS results (top panels) show the incoherent nature 
of the flows — it is this that leads to the suppression of the jets. For this case, the magnetic energy 
is in approximate equipartition with the kinetic energy, as shown in Table 1. Once again DNS and 
DSS show a remarkable agreement for this case; see Fig. [6j A comparison of the mean toroidal 
component of the magnetic field is shown in Fig. [7j Here the situation is reversed from the previous 



9 We note that sufficient averaging must be employed in order to obtain meaningful averages in both DSS and 
DNS. This is easy to achieve for DSS but is problematic for DNS. 
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weaker field case. The magnetic field here is too strong to be expelled by the eddies and the jet 
never forms at low latitudes. Therefore the field is confined to low latitudes and the (weaker) jets 
can only be found at high latitudes where the imposed field is weaker. For this strength of imposed 
field the kinetic energy is reduced (see the values of Rm in Table 1) and the magnetic energy 
comes into equipartition with the kinetic energy. Clearly the transport of angular momentum by 
the Reynolds and Maxwell stresses and of magnetic flux by the turbulent advection has acted in a 
very different manner here. The discussion of these processes and their description via the second 
cumulants is postponed to a subsequent paper. 



5.2. Case L m i n — 8, M m i n — 1, and L max — 10 

We now consider the effect of reducing the minimum stochastic forcing wavevector in the 
azimuthal direction, M m j n , down to wavenumber 1. This brings stochastic effects to larger scales 
and so presents a more robust challenge for DSS. A comparison of the zonal mean relative vorticity 
and zonal velocity as calculated by DNS and DSS (CE2) is shown in Fig. [8] for the hydrodynamic 
problem. In contrast to the previous case, there are two prograde jets at high latitudes, and one 
equatorial retrograde jet; again the total angular momentum is close to zero. Now, however, the 
jets are seen to wander significantly in latitude in DNS owing to the continual random forcing at 
large zonal scales. Once established in DSS, however, they remain fixed in place. As a consequence, 
the time-averaged zonal means are reduced in magnitude in DNS compared with DSS, as made 
apparent in the quantitative plot of Fig. [9j 

Imposing a relatively weak toroidal field by setting Bq = 0.1 again has little effect on the zonal 
mean velocity as shown in Figure [9j However at Bq = 1 the jets are largely eliminated by the 
fluctuations in the magnetic field, both in DNS and in DSS. Somewhat larger jet remnants remain 
at high latitudes, where the imposed toroidal field is weakest, and again stronger jets are found in 



DSS than in DNS (see Fig. 10) 



The toroidal field is tightly confined to latitudes less than roughly 60°, as Fig. [TT] depicts, 
likely owing to flux expulsion of the field. Again DSS does a reasonably good job of reproducing 
the mean field. 



6. Discussion 

In this paper we have introduced the concept of Direct Statistical Simulation (DSS) for astro- 
physical fluid dynamics. We have compared the results of Direct Numerical Simulation (DNS) and 
DSS for the problem of two-dimensional hydrodynamics and magnetohydrodynamics on a spherical 
surface. Although the set-up of the model is relatively simple, the ensuing dynamics is not. In 
the hydrodynamic case, non-trivial interactions at moderate scales drive inhomogeneous large-scale 
zonal flows (jets/winds). With a weak imposed field the jets remain largely unaffected and the 
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magnetic fields are expelled to higher latitudes. With a stronger imposed toroidal magnetic field 
these winds are suppressed except at high latitudes where the imposed field is weak. 

We find that even the simplest formalism of DSS, based upon the truncation of the cumulant 
hierarchy at second order, is capable of reproducing the driving and suppression of the zonal flows 
and the flux expulsion of the magnetic fields by the inhomogeneous jets. Because the method 
includes interactions that are non-local in space it is very well suited to such inhomogeneous prob- 
lems thatypically arise in astrophysics. Such a truncation is equivalent to keeping the mean/eddy 
interactions in the eddy equations and the eddy/eddy interactions in the mean equations, whilst 
suppressing the eddy/eddy interactions in the eddy equations. Thus it is the Reynolds and Maxwell 
stresses that respectively drive and suppress the jet, not an inverse cascade as is frequently assumed. 

The DSS scheme is more numerically efficient than the corresponding DNS. We believe that the 
results presented here are an encouraging beginning for the concept of DSS in astrophysical fluid 
dynamics. It is important though to determine the range of validity of such a procedure. Clearly the 
method is designed to work best when the dynamics leads to the generation of substantial statistical 
means (e.g. mean flows or magnetic fields) or involves the interactions of prescribed (usually 
inhomogeneous) mean quantities with smaller scale turbulence. The method is inefficient when 
the turbulence is dominated completely by small-scales and is largely homogeneous — for example 
homogeneous, isotropic hydrodynamic turbulence, or the small-scale dynamo problems (see e.g. 
Tobias et al.|2011 ). We do believe, however, that many cases of astrophysical interest do fall into the 



category where DSS techniques may prove useful. Examples currently under consideration include 
the interaction of mean magnetic fields and shear flows either on a spherical surface (leading to 
joint instability (see e.g. Cally et al.|20 03)) or in a cylindrical domain (leading to magnetorotational 



instability), the instability and mixing of large-scale shear flows in the presence of a magnetic field, 



and the driving of zonal flows via convection in a tilted cylindrical annulus (see e.g. Brummell & 
Hart||1993 ). It will be interesting to determine how well the techniques described in this paper fare 
for these problems, and we predict varying degrees of success. We also stress that even utilising 
DSS will not allow the calculation of statistics at astrophysically realistic values. However it is to 
be hoped that, whereas the dynamics may be extremely sensitive to the parameters, for a range 
of problems the statistics may prove less so. We believe that some of these problems may require 
the inclusion of higher order cumulants in the scheme and are currently engaged in determining 
efficient numerical procedures for their integration. 

Clearly in the longer term, if these techniques prove useful for the simpler problems described 
above, it will be of interest to apply them to more computationally intensive problems in astrophys- 
ical fluids. These include the driving of zonal flows in planets, the mixing of angular momentum 
and abundances in stellar interiors and the transport by turbulence in accretion disks. 
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Case 


Method 


B 


Ro 


Rm 


(B 2 )/(U 2 ) 


1 


DNS 





0.0351 


4415 


0.0 


1 


DSS(CE2) 





0.0352 


4431 


0.0 


1 


DNS 


0.1 


0.0331 


4195 


0.04 


1 


DSS(CE2) 


0.1 


0.0323 


4064 


0.04 


1 


DNS 


0.5 


0.015 


1865 


1.29 


1 


DSS(CE2) 


0.5 


0.017 


2144 


1.02 


2 


DNS 





0.0542 


6812 


0.0 


2 


DSS(CE2) 





0.0548 


6885 


0.0 


2 


DNS 


1.0 


0.0237 


2980 


1.07 


2 


DSS(CE2) 


1.0 


0.0396 


4980 


1.25 



Table 1: Non-dimensional numbers as calculated a posteriori. Case 1 refers to stochastic forcing 
with M m i n = 8 whilst Case 2 refers to M m j n = 1. (For the case of zero magnetic field, Rm reduces 
to a non-dimensional measure of the kinetic energy of the flow.) 




Fig. 1. — DNS calculation of the instantaneous relative vorticity (left panels) and zonal velocity 
fields (right). Top: Pure hydrodynamic model (Bq = 0) that develops a westerly (prograde) jet 
along the equator. Bottom: Imposed toroidal magnetic field parameter Bq = 0.5. The calculation 
is done on a spherical geodesic grid with D = 163, 842 cells with stochastic forcing over spherical 
wavevectors 8 < £, m < 12. See Sec. 4.2 for the values of the other parameters. 
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Fig. 2. — Timelines of the zonal mean relative vorticity (left) and zonal velocity (right) as calculated 
for pure hydrodynamic problem with no imposed toroidal field (Bq = 0). Top: DNS on a spherical 
geodesic grid with D = 163, 842 cells. Bottom: DSS (CE2) with L = 100 (bottom). A running time 
average commences at the midpoint in time (vertical line) at which point statistics are accumulated. 
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Fig. 3. — Comparison of the mean zonal velocity as calculated in DNS and DSS (CE2) in the pure 
hydrodynamic problem with no imposed magnetic field (Bq = 0). Convergence with increasing 
resolution is evident both for DNS and for DSS (CE2). The prograde jet is reproduced well by 
CE2. Due to the finite time interval of accumulating statistics (1000 days), statistics obtained from 
DNS are not perfectly symmetric about the equator. 
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Fig. 4. — Timelines of the zonal mean zonal velocity for an imposed toroidal magnetic field with 
B = 0.1. Top: DNS with D = 40,962 cells. Bottom: DSS (CE2) with L = 100. 




Fig. 5. — Same as Fig. [2] except with Bq = 0.5. 
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Fig. 6. — Comparison of the mean zonal velocity as calculated in DNS and DSS (CE2) for imposed 
toroidal fields of Bq = 0, 0.1 and 0.5. 




Fig. 7. — Comparison of the mean toroidal magnetic field as calculated in DNS and DSS (CE2) for 
imposed toroidal fields of Bq = 0.1 and 0.5. 
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Fig. 8. — Timelines of the zonal mean relative vorticity (left) and zonal velocity (right) as calculated 
for pure hydrodynamic problem with no imposed toroidal field (Bq = 0). Top: DNS on a spherical 
geodesic grid with D = 163,842 cells. Bottom: DSS (CE2) with L = 100 (bottom). 




Fig. 9. — Comparison of the mean zonal velocity as calculated in DNS and DSS (CE2) for imposed 
toroidal fields of Bq = and 1.0. 
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Fig. 10. — Same as Fig. [8] except for B = 1 and with DNS run on a lower resolution spherical 
geodesic grid with D = 40, 962 cells. 




Fig. 11. — Toroidal component of the magnetic field for Bq = 1.0. Top: Instantaneous field found 
by DNS on a spherical geodesic grid with D = 40, 962 cells. Bottom: Zonal mean field found by 
DSS with L = 50. 



